function [sgb,equalorder, ...
          Ugridreduced, sgreduced,...
          id_Ugridreduced,...
          s_id_gridreduced, sU_id_gridreduced,...
          n_U_sets, Tcoord, indices_pairs,...
          numeqconstraints,coordrowseq, ...
          d3, e5, f5, g3, h5, i5, l3, m5, n5,...
          fillAeq, beq]=useful_anyparam3_ind(u01grid, u11grid, u21grid,u02grid, u12grid, u22grid, sg,PYX_cond)     
%% U sets
n_U_sets=3; %number of U-sets (number of dimensions of the distribution of the taste shock differences)
Ugrid=cell(1,n_U_sets); %cell 1xn_U_sets
Ugrid{1}=[u01grid-u11grid -u01grid+u11grid u01grid-u21grid -u01grid+u21grid u21grid-u11grid -u21grid+u11grid ...
         u02grid-u12grid -u02grid+u12grid u02grid-u22grid -u02grid+u22grid u22grid-u12grid -u22grid+u12grid ...
         zeros(sg,1) Inf*ones(sg,1) -Inf*ones(sg,1)]; 
Ugrid{2}=[u01grid-u21grid -u01grid+u21grid u01grid-u11grid -u01grid+u11grid u21grid-u11grid -u21grid+u11grid ...
         u02grid-u22grid -u02grid+u22grid u02grid-u12grid -u02grid+u12grid u22grid-u12grid -u22grid+u12grid ...
         zeros(sg,1) Inf*ones(sg,1) -Inf*ones(sg,1)]; 
Ugrid{3}=[Inf*ones(sg,1) -Inf*ones(sg,1) ...
         u21grid-u11grid -u21grid+u11grid u01grid-u11grid -u01grid+u11grid u01grid-u21grid -u01grid+u21grid ...
         u22grid-u12grid -u22grid+u12grid u02grid-u12grid -u02grid+u12grid u02grid-u22grid -u02grid+u22grid ...
         zeros(sg,1)]; 
sgb=0;     

%% Reduce number of grid points to check by applying our partitioning arguments
%equalorder (sgx1) assigns the same index to parameters vectors such that Ugrid{1}, Ugrid{2}, ..., Ugrid{n_U_sets} have the same order

%%%%%%%% Criterion 1 %%%%%%%%
d=cell(n_U_sets,1);
index=cell(n_U_sets,1);
for j=1:n_U_sets 
    [sortedUjgrid,index{j}] = sort(Ugrid{j},2); 
    %sortedUjgrid gives Ugrid{j} with each row sorted from smallest to largest; it is sgxn_columns_Ujgrid
    %for each row i of sortedUjgrid and for k=1,2,...,n_columns_Ujgrid, index{j}(i,k) gives the column position in Ugrid{j}(i,:) of sortedUjgrid(i,k); it is sgxn_columns_Ujgrid
    d{j} = diff(sortedUjgrid,[],2) == 0; %sg x (n_columns_Ujgrid-1)
    %Suppose that Ugrid{j} has 3 columns. Then, for each row i of sortedUjgrid, d{j} will be: 
    %[1 1] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[1 0] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
    %[0 1] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[0 0] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
end
%%%%%%%% Criterion 2 %%%%%%%%
n_columns_Ujgrid=15;
Ugrid_sum=zeros(sg,n_columns_Ujgrid+n_columns_Ujgrid^2);
for g=1:sg
Ugrid_sum(g,:)=[Ugrid{1}(g,:) reshape((Ugrid{2}(g,:)+Ugrid{3}(g,:).')',[],1)'];
end

[sortedUgrid_sum,index_sum] = sort(Ugrid_sum,2); 
d_sum = diff(sortedUgrid_sum,[],2) == 0; %sg x (n_columns_Ujgrid*n_U_sets-1)

%%%%%%%% Combine %%%%%%%%
[~,~,equalorder] = unique([index_sum d_sum [index{:}] [d{:}] ],'rows', 'stable'); %sgx1 

%Reduced grid of parameter vectors and U sets
[~,index_final,~]=unique(equalorder, 'stable'); %one row from equalorder per ordering group together with its index
Ugridreduced=cell(1,n_U_sets);
for j=1:n_U_sets 
    Ugridreduced{j}=Ugrid{j}(index_final,:); %sgreduced x n_columns_Ujgrid
end
sgreduced=size(Ugridreduced{j},1); %size of the reduced grid of parameter values



%% Construct indices to keep track of equal elements inside Ugridreduced{j} for j=1,...,n_U_sets
id_Ugridreduced=cell(n_U_sets,1);
for j=1:n_U_sets 
    id_Ugridreduced{j} = NaN(size(Ugridreduced{j})); % %sgreduced x n_columns_Ujgrid
    for k = 1:sgreduced
        [~, ~, id_Ugridreduced{j}(k,:)] = unique(Ugridreduced{j}(k,:), 'stable'); 
    end
end

s_id_gridreduced=zeros(sgreduced, n_U_sets); %useful to set number of unknowns of the linear programming
for j=1:n_U_sets
    s_id_gridreduced(:,j)=max(id_Ugridreduced{j},[],2); %number of distinct elements of Ugridreduced{j}
end
sU_id_gridreduced=prod(s_id_gridreduced,2); %sgreducedx1 reporting the cardinality of the Cartesian product of the U sets
                                            %sU_id_gridreduced(g) is the number of unknowns of the LP for each g-th parameter value
                                            
s_id_gridreduced_2=ones(1, n_U_sets)*n_columns_Ujgrid; %number of elements of Ugridreduced{j} (count as separate even equal elements)
sU_id_gridreduced_2=prod(s_id_gridreduced_2); %cardinality of the Cartesian product of the U sets     
%% Invariant components of the linear programming problem 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Inequality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vectors = arrayfun(@(x) {1:x}, s_id_gridreduced_2); 
n = numel(vectors);
Tcoord_temp = cell(1,n);
[Tcoord_temp{:}] = ndgrid(vectors{:}); 
Tcoord_temp = cat(n+1, Tcoord_temp{:});
Tcoord = reshape(Tcoord_temp,[],n); %matrix of dimension sU_id_gridreduced_2xn_U_sets reporting the column coordinates of all the possible n_U_sets-tuples from the U-sets
%For example if n_U_sets=3 and s_id_gridreduced_2 is [2 3 1], then 
%[ca, cb, cc] = ndgrid(1:2, 1:3, 1:1);
%Tcoord = [ca(:), cb(:), cc(:)];   
%Tcoord= [1     1     1
%         2     1     1
%         1     2     1
%         2     2     1
%         1     3     1
%         2     3     1];
indices_pairs=pairIndices(sU_id_gridreduced_2); %row indices of all possible unordered pairs of rows from the matrix obtained by taking the Cartesian product of the U sets
                                         %[sU_id_gridreduced_2*(sU_id_gridreduced_2-1)/2]x[2]   
%Continuing the example above
%indices_pairs=[1 1; 1 2; ...; 1 27; 2 3; ...; 2 27; ... 26 27]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Equality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

numeqconstraints=6+...
                 (15+3*2*sgb+15+3*2*sgb-1)*(15+3*2*sgb-1)+(15+3*2*sgb)^2+...
                 1+...
                 (6+1+sgb*3)*3+...
                 (2*13)+(2*2*sgb)*3;  
             
numzeros=(15+3*2*sgb+15+3*2*sgb-1)*(15+3*2*sgb-1)+(15+3*2*sgb)^2;

coordrowseq=[1 1 1 ... %observational equivalence
             2 2   ...
             3     ...
             4 4 4 ...
             5 5  ...
             6     ...
             7:1:6+numzeros ... %zeros
             6+numzeros+1 ... %ones
             kron(6+numzeros+1+1:1:6+numzeros+1+(6+sgb*3)*3, ones(1,2)) ... %symmetry
             6+numzeros+1+(6+sgb*3)*3+1:1:6+numzeros+1+(6+sgb*3)*3+3 ...
             kron((6+numzeros+1+(6+sgb*3)*3+3+1:1:6+numzeros+1+(6+sgb*3)*3+3+(2*13)+(2*2*sgb)*3),ones(1,2))]; %identical  

         
%used to construct coordcolumnseq
d1=kron((16:1:15+2*sgb), ones(1,2));
d2=14*ones(1,2*sgb*2);
d3=[d1;d2];

e1=15+2*sgb+1:1:15+2*sgb+2*sgb;
e2=14*ones(1,2*sgb);
e3=[e1;e2];
e3=e3(:).';
e4=14*ones(1,2*sgb+2*sgb);
e5=[e4; e3];

f1=ones(1,2*sgb);
f2=15+2*sgb+1:1:15+2*sgb+2*sgb;
f3=[f1;f2];
f3=f3(:).';
f4=ones(1,2*sgb+2*sgb);
f5=[f4; f3];

g1=kron(15+2*sgb+1:1:15+2*sgb+2*sgb, ones(1,2));
g2=14*ones(1,2*sgb+2*sgb);
g3=[g1;g2];

h1=16:1:15+2*sgb;
h2=14*ones(1,2*sgb);
h3=[h1;h2];
h3=h3(:).';
h4=14*ones(1,2*sgb+2*sgb);
h5=[h4; h3];

i1=ones(1,2*sgb);
i2=15+4*sgb+1:1:15+4*sgb+2*sgb;
i3=[i1;i2];
i3=i3(:).';
i4=ones(1,2*sgb+2*sgb);
i5=[i4;i3];

l1=kron((15+4*sgb+1:1:15+4*sgb+2*sgb), ones(1,2));
l2=14*ones(1,2*sgb+2*sgb);
l3=[l1;l2];

m1=15+4*sgb+1:1:15+4*sgb+2*sgb;
m2=14*ones(1,2*sgb);
m3=[m1;m2];
m3=m3(:).';
m4=14*ones(1,2*sgb+2*sgb);
m5=[m4;m3];

n1=ones(1,2*sgb);
n2=16:1:15+2*sgb;
n3=[n1;n2];
n3=n3(:).';
n4=ones(1,2*sgb+2*sgb);
n5=[n4;n3];


                         
fillAeq=[1 -1 -1 ... %observational equivalence
         1 -1    ...
         1       ...
         1 -1 -1 ... 
         1 -1    ...
         1       ...
         ones(1,numzeros) ... %zeros
         1 ...%ones 
         repmat([1 1],1, (6+sgb*3)*3) ...%symmetry
         ones(1,3)...
         repmat([1 -1], 1, (2*13)+(2*2*sgb)*3)]; %identical  
            
beq=[ PYX_cond(1)-1;... %observational equivalence
      PYX_cond(2);...
      PYX_cond(3);...
      PYX_cond(4)-1;...
      PYX_cond(5);...
      PYX_cond(6);...
      zeros(numzeros,1);... %zeros
      1;... %ones       
      ones((6+sgb*3)*3,1); ... %symmetry
      1/2*ones(3,1)
      zeros((2*13)+(2*2*sgb)*3,1)]; %identical  %[#equalityconstraints]x1       
end